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Abstract. To investigate the stability properties of polar disks we performed two-dimensional hydrodynamical simulations for 
flat polytropic gaseous self-gravitating disks which were perturbed by a central SO-like component. Our disk was constructed 
to resemble that of the proto-typical galaxy NGC 4650A. This central perturbation induces initially a stationary two-armed 
tightly- wound leading spiral in the polar disk. For a hot disk (Toomre parameter Q > 1.7), the structure does not change over 
the simulation time of 4.5 Gyr. In case of colder disks, the self-gravity of the spiral becomes dominant, it decouples from the 
central perturbation and grows, until reaching a saturation stage in which an open trailing spiral is formed, rather similar to that 
observed in NGC 4650A. The timescale for developing non-linear structures is 1-2 Gyr; saturation is reached within 2-3 Gyr. 
The main parameter controlling the structure formation is the Toomre parameter. The results are surprisingly insensitive to the 
properties of the central component. If the polar disk is much less massive than that in NGC 4650A, it forms a weaker tightly- 
wound spiral, similar to that seen in dust absorption in the dust disk of NGC 2787. Our results are derived for a polytropic 
equation of state, but appear to be generic as the adiabatic exponent is varied between y = I (isothermal) and -y = 2 (very stiff). 
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1. Introduction 

Polar ring galaxies (PRGs) are characterized by an early-type 
_ or elliptical central component and a gas-rich ring or disk 
■ which is nearly polar with respect to the central galaxy. The 
' spatial orientation and the kinematics of the main components 
' makes them "multi-spin systems" (Rubin 1994. 1 with mutually 
' orthogonal spin vectors. In the polar-ring catalog (Whitmore 
et al.ll990l about 100 PRGs or PRG candidates are given. The 
polar structures can be either naiTOW rings, or radially-extended 
disks. Both morphologies are about equally frequent. 

A prototype of a PRG with a radially-extended disk is NGC 
4650A. It is characterized by a central gas-poor SO compo- 
nent and an exponential polar disk much larger than the cen- 
tral galaxy. At an adopted distance of 35 Mpc, the central 
component has a luminosity Lj ^ 3 ■ lO' L0, coiTespond- 
ing to a mass of 5 ■ 10'' M0 in stars, while for the polar disk 
L/ =s 1.2 ■ 10'' Lq, so the stellar mass 3 ■ lO'' Mq (Gallagher et 
al. I2OO2I 1. The polar disk shows roughly an exponential dis- 
tribution of stellar light (lodice et al. I2002I I. It is gas-rich, 
with an Hl-mass of about 7-10^ Mq at the adopted distance 
(Arnaboldi et al. 1 997j . and about 10% as much mass in molec- 
ular form (Watson et al. ll994> . The HI maps show gas rotating 
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at speed 120kms ' outtoaradius of25kpc, which corresponds 
to a dynamical mass of 8 ■ 10'" Mq. The polar disk is blue, with 
about 1 M0 yr ' of ongoing star formation, whereas in the SO 
component star formation appears to have ceased 3-5 Gyrs 
ago (Gallagher et al. 2002 1. 

Arnaboldi et al. pointed out the presence of two open spi- 
ral arms in the polar disk, with approximate mirror symmetry 
about the central SO, in both optical and near-infrared light. 
These have a spatial extent of about 20 kpc in diameter (see 
their Fig. 10 or our Fig. ll3lbelow). The arms are also traced 
by a few blue knots, suggesting that they represent real density 
enhancements where stars have formed, and not simply a warp 
in the polar disk (Gallagher et al. 2002). 

The polar disk in NGC 2787, which is classified SBO/a with 
a LINER nucleus, exhibits a very different morphology, with 
a tightly-wound spiral in the inner 7 arcsec or 200 pc of the 
polar disk. The HST/WFPC2 image (see Figure 13 of Erwin 
& Sparke 2003jor Figure 4 of Sil'chenko & Afansiev l2004t 
shows the spiral in dust absorption SB of the nucleus where 
the polar disk lies in front of the galaxy. The central galaxy is 
about as luminous as NGC 4650A, with a; 2 ■ 10^ L0 at a 
distance of 7.5 Mpc, and probably has a similar mass. But the 
polar structure is much less massive, with only 3-10'* M0 of 
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HI and H2 (Welch & Sagel2003'l, and most of this is in a ring 
at radius ~ 6.4 arcmin or 14 kpc. 

Though polar ring galaxies (PRGs) are not very frequent 
(about 0.5% of all galaxies (Whitmore et al. |1990|f), they at- 
tracted the interest of astronomers mainly for two reasons. 
First, the fact that they are multi-spin systems makes them 
prime targets to study rotation curves in orthogonal directions. 
From that, the flattening of dark haloes might be deduced 
(e.g. Whitmore et al. 1987 Sackett & Sparke Sackett 
et al. [T994]) . Second, their formation mechanism is still un- 
der debate though it is assumed that galaxy interactions play 
a key role. Basically two formation concepts have been dis- 
cussed in the past. One scenario is a tidal accretion model, in 
which mass is transferred in an encounter between two galaxies 
(e.g. Reshetnikov & Sotnikova 1997 1. The second assumes a 
head-on merger between two orthogonal spiral galaxies (Bekki 
1199 8 1. Bournaud & Combes ( 2003 1 performed numerical simu- 
lations to compare the probability of both scenarios, favouring 
the accretion scenario. However, a lot of assumptions which 
are observationally hard to determine enter the estimate of the 
probabilities for both scenarios. 

In order to understand the formation and evolution of polar 
disks, their ages and, thus, their stability properties are of spe- 
cial interest. From their regular morphology, many appear to be 
stable on at least several dynamical timescales (see e.g. Sparke 
120041 for a review). If so, differential precession effects have to 
be suppressed effectively. Steiman-Cameron & Durisen ( 1982) 
suggested that a non-spherical halo might be a stabilizing fac- 
tor. Sparke {1986 ) and Arnaboldi & Sparke J1994> emphasized 
that self-gravity may stabilize polar disks. Especially, in mas- 
sive polar disks like in NGC 4650A self-gravity is expected to 
be important. 

In this paper, we want to investigate the stability proper- 
ties of self-gravitating polar disks which are perturbed by a 
central galactic component. Our numerical analysis is based 
on two-dimensional simulations for flat disks. Such calcula- 
tions have been widely used for modelling stellar and gaseous 
disks (e.g. Englmaier & Shlosman 2000 Orlova et al. 2002 
Masset & Bureau 12003b . for simulating star-forming disks 
(e.g. Korchagin & Theis 119991 or dusty nuclear disks (Theis 
& Orlova 2004). This numerical approach allows not only to 
study the evolution during the linear stage, but also to follow 
the system deeply into the non-linear regime. From that we get 
e.g. the timescales of the different evolutionary stages or the 
spatial pattern and kinematics of the system. 

The numerical code is described in Sect. |2] In Sect. |3l the 
initial model is presented. The numerical simulations are given 
in Sect.|3and|5] Finally, we discuss and summarize our results 
in Sect.|6l 

2. Hydrodynamic code 

2.1. Numerical Method 

In our simulations we solve the two-dimensional hydrodynam- 
ical equations for a flat disk in cylindrical coordinates (with 
z - 0). This represents a polar ring or disk around a galaxy that 
is flattened in the x-z plane, as shown in Figure^ In some im- 



ages cartesian coordinates are shown: then the positive x-axis 
corresponds to = and tf increases counter-clockwise. We 
solve the 
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2. the Euler equation for the radial velocity u 
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The gravitational potential <1) = Ojisk + 'I'hb + ^bai consists of 
three components: the polar disk, a halo/bulge, and a contri- 
bution from the flattened SO galaxy. The gravitational poten- 
tial Odisk of the flat polar disk is given by (Binney & Tremaine 
IW7I hereafter BT87, Sect. 2.8) 
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This sum is calculated using a FFT technique applying the 
2D Fourier convolution theorem for polar coordinates. External 
potentials describe the contribution of a rigid halo/bulge. The 
latter can be done in two ways: either a halo/bulge mass dis- 
tribution is given from which the equilibrium rotation curve 
is derived (considering also the pressure forces and the self- 
gravity of the disk) or the rotation curve vdR) is given and the 
halo/bulge mass distribution is derived. In the latter case (which 
we used here), the contribution d(t>hb/dR, i.e. the enclosed mass 
of the halo/bulge component, is calculated from 



v^JR)^R-^((t> + <:>bb). 
dR 



(5) 



Though formally, any profile of Ohb might be a solution of Eq. 
(|5}, we only accept physically meaningful, radially increasing 
cumulative mass distributions. More details of our setup of the 
rotation curve can be found in Sect. 13.21 

In the X -y plane of the polar material, the oblate SO galaxy 
is represented by the bar-like potential of a stationary Ferrers' 
bar (BT87, Sect. 2.3). In order to avoid an initial "kick" on 
the disk, the mass Mbai of the bar-like perturbation is slowly 
enhanced starting from M\,^(t = 0) = 0. During this adiabatic 
turn-on process the mass Mhb,<a of the halo/bulge inside the 
semi-major axis a is reduced by keeping the sum of the masses 

Mbar(f) + Mhb,<a(0 = COnSt. 

For the pressure terms we applied a polytropic equation of 
state 



P^C-IJ 



(6) 



with an adiabatic exponent of y = 5/3. This allows to close the 
set of hydrodynamical equations with the equations for the first 
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order moments, i.e. the Euler equations (|5Jl and Q. The hy- 
drodynamic equations Q - (|3} are discretized on a logarithmic 
Eulerian grid with a constant logarithmic grid spacing. The dif- 
ferent terms in the equations Q - Q are taken into account by 
applying an operator splitting technique like e.g. in the ZEUS 
program (Stone & Norman 1992 1. Advection is performed by 
the second-order van-Leer scheme. The timestep is determined 
according to the usual Courant-Friedrichs-Levy criterion. For 
our simulations we use reflecting boundary conditions in radial 
direction (and periodic in azimuthal direction). Our simulations 
are mainly done on a grid with 270x270 cells. 

Most of the simulations are performed without any ex- 
plicit artificial viscosity (AV). Test calculations with differ- 
ent amounts of AV did not significantly differ from the refer- 
ence simulation without AV. During the phase of linear growth 
the global amplitudes show almost no quantitative difference. 
Later, in the saturation stage, shocks become more important 
and AV gives different results. However, the saturation levels 
just scatter around the models without AV. Because of that and 
because we are mainly interested in the growth time for an 
instability, we did not apply AV in most of the simulations, 
though it might improve the treatment of shocks in the late sat- 
uration stage. 

2.2. Units 

The units are chosen to be 10^ Mq and 1 kpc for the mass and 
length scale, respectively. The gravitational constant was set 
to G = 1. This results in a velocity unit of 65.6 kms ' and a 
unit for the angular speed of 65.6 kms ' kpc The time unit 
is 14.9 Myr. These units are used throughout the paper, unless 
other units are given explicitly. 

2.3. Dimensions, Times and Accuracy 

Motivated by the polar disk of NGC 4650A we adopted an in- 
ner edge of 2 kpc for most of our simulations. The outer edge 
is set to 50 kpc to prevent perturbations generated at the outer 
edge from affecting the solutions in the "central" area of the 
polar disk (10 kpc). When using a total of 270 radial grid cells, 
our logarithmically equidistant grid has a radial cell size of 24 
pc at the inner edge degrading to 0.6 kpc at the outer boundary. 

A typical simulation run needs of the order of several 
10^ to 10^ timesteps, until the end of the calculation at 
t = 300 ~ 4.47 Gyr is reached. A typical timestep is of the 
order of Af ~ 1 ... 6 ■ 10"^ ~ 1 .5 ... 9 ■ 10"* yr. The simulations 
were performed on a NEC-SX5 vector computer The code 
is highly vectorized, reaching up to 50% of the peak perfor- 
mance of a single vector processor (sustained performance: 
1.6-2 GFlop). The most time-consuming part of the simula- 
tions is the calculation of the disk's self-gravity. Numerically 
this is done with a fast Fourier transform from the MathKaisan 
library on the NEC-SX5 computer which is especially fast for 
this number of grid cells. The usual number in that resolution 
regime is 256 = 2^ cells per dimension. However, the speed 
of the FFT becomes rather slow for this array size due to bank 
conflicts in the memory access. 
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Fig. 1. Rotation curve vv and contributions of the different com- 
ponents (v'hb: halo&bulge, Vdisk: polar disk and vp: pressure) for 
the unperturbed model A. Note that the contributions are added 
up quadratically for v^. 

For the simulation of the unperturbed disk (model A) the 
energy was conserved to better than 5 ■ 10"^ (relative deviation), 
angular momentum better than 2-10"^ and the mass even within 
the machine accuracy of double-precision numbers. The accu- 
racy of the code decreases when the system becomes non-linear 
and the saturation regime is reached. In that case the accuracy 
is of percent level for both, energy and angular momentum. 

3. Initial model 

In order to mimic real polar disks to a reasonable degree the 
parameters for our initial model were chosen to match closely 
the observations of the prototype polar disk NGC4650A 
(Gallagher et al. '2002"(. However, we do not intend to have a 
near-to-exact model for NGC4650A, but a "prototypical" po- 
lar disk. Thus we preferred to describe the rotation curve or 
the mass distribution by "simple" standard formulas instead of 
matching exactly small- or intermediate-scale variations. We 
also do not consider the multi-phase nature of the ISM or (en- 
ergy feedback from) star formation in our calculations. In detail 
the model is characterized by the mass distribution, the kine- 
matics of the disk and the appUed equation of state. 

3.1. Poiardisk 

In our standard model, a mass of Mdisk = 1-2 ■ 10'° M0 is 
distributed exponentially with a scale length of 4 kpc. The 
"punched" disk has a sharp inner boundary at a radius of 2 kpc. 
The half-mass radius of the disk is about 7.3 kpc, and 90% 
of its mass lies within the 'central' region at < 16 kpc. 
Additionally, the surface density is slightly perturbed result- 
ing in a random deviation from the mean surface density of an 
amplitude of 10"^. 

3.2. Rotation curve 

In principle our code allows us to adopt any prescribed rota- 
tion curve by adjusting the halo/bulge mass distribution accord- 
ingly. In order to determine the required halo/bulge potentials 
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the gravity of the disk as well as the pressure gradient terms 
are taken into account. They are subtracted from the rotation 
curve in order to get the remaining halo contribution. Though 
this scheme allows always for a formal solution, only physical 
cases are accepted by requiring that the enclosed halo/bulge 
mass is monotonically increasing with radius. The halo/bulge 
potential is non-rotating. The bar-like SO component is treated 
as stationary except for a short, merely technical "adiabatic" 
start phase (which is described below). 

For the models in this paper we parametrized the rotation 
curve by 



VciR) 



(7) 



(10) 



This guarantees a rigid rotation in the inner parts and a flat rota- 
tion curve outside. The parameter n defines the smoothness of 
the transition between these two zones. For our simulations we 
mainly used a sharp transition near Rn-^t = 2 kpc by choosing 
n = 10. Hence, the rotation curve is basically flat all over the 
polar disk (Fig.Q. 

The initial azimuthal velocity of the gaseous phase is calcu- 
lated by the Jeans' equation for the radial velocity component 
which reads in cylindrical coordinates 

du du V du ^ dP d , ^ ^ ^ 

— -I- M — + — ■ (0 + Ohb) • (8) 

dt dR Rd4> R 2 dR dR 

In case of the initial equilibrium, the radial velocity u vanishes. 
If we define the pressure contribution vp to the rotation curve 
as 

2 _ R dP 

= "z ■ M ' 

then the azimuthal velocity v is given by 
d R dP . 

^'^m^'^^ '^''^ ^^"dR^ '^'^'^ ~ ''^^^ ■ 

The rotation curve is derived from the gravitational potential 
according to Eq. (|5}- 

Decomposing disk and halo contributions, we find that our 
disk is submaximal by a factor of 3.6. This leaves suflicient 
mass for an additional central component. It should be noted 
that close to the inner edge the gravitational pull of the disk 
becomes negative, i.e. points outwards. This "repulsive" force 
has to be compensated for by a halo contribution corresponding 
to a rotational velocity exceeding the rotation curve. 

We choose the rotation speed Voa at large radii to be 
150 kms The rotation period at the inner edge is 5.9 time 
units (^ 90 Myr) increasing to ~ 300 Myr at the half-mass ra- 
dius of the disk and reaching ~ 2 Gyr at its outer edge at 50 kpc. 
The maximum of Q.{R) - k{R)/2 is important for the existence 
of inner Lindblad resonances: here k is the epicyclic frequency 
and Q = v^/R the circular speed. This maximum is reached 
at 2.55 kpc where D.-k/2^ 0.236 or = 15.4 kms"' kpc"'. In 
case of pure stellar disks, a two-armed spiral pattern rotating at 
the angular rate Q.p can propagate only where \Q.p - f2| < k/2, 
so any such spiral with Qp < 0.236 cannot propagate from the 
outer disk to the inner boundary. In our study we deal, how- 
ever, with gaseous disks, in which density waves can propa- 
gate as sound waves through the Lindblad resonances, with 




Fig. 2. Initial radial profile of the Toomre parameter Q for the 
unperturbed model A and the reference model B. The minimum 
2 of 1.4 is reached at about 6 kpc. 



little damping (for a brief summary see Bertin (20001, here- 
after GBOO, Chap. 16). Polar rings are extremely gas-rich; of- 
ten more of their baryonic mass is in cool gas than in stars (e.g. 
Sparke l2004> . So this is likely to be an appropriate approxima- 
tion. 

3.3. Toomre parameter 

The local stability of a disk with respect to axisymmetric per- 
turbations is controlled by the Toomre parameter 



(9) 



ttGI. 



(11) 



Q compares the stabilizing factors - velocity dispersion ctr 
(sound speed) and radial oscillation frequency (or epicyclic fre- 
quency k) - with the destabilizing role of the self-gravity of a 
mass element with surface density S (e.g. BT87 Sect. 6.2). A 
value of Q less than unity means instability to axisymmetric 
perturbations. In gaseous disks the most unstable wavelength 
A„i in the radial direction is given by 



with the critical wavelength A^nt defined as 



Arrtt = 



(12) 



(13) 



From numerical experiments it is known that disks become 
rather stable even to non-axisymmetric modes for minimum Q- 
values above 1 .4. In case of flat rotation curves one can derive 
from linear stability analysis that Q > V3 ~ 1.73 guaran- 
tees that a fluid disk is stable with respect to all perturbation 
modes (see e.g. Bertin et al. ( l_989j, Polyachenko et al. ( 1997 i 
or Section 15.2 of GBOO). This is also in good agreement with 
multi-component disks studied by Orlova et al. (2002'l. In our 
reference model the constant C of Eq. (|6} is adjusted in such 
a way that the initial minimum value of 2 is 1 .4. This value is 
reached at R x 6 kpc. Fig.|3shows the radial run of Q for this 
model. 
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Table 1. Properties of the numerical models 



time = 300.0 



model 


MU Mq) 


Qmm 


comment 


A 




1.4 


no central component 


B 


5 


10^ 


1.4 


reference model 
(with central component) 


MBl 


1 


10'^ 


1.4 


small SO component 


MB2 


1 • 


10'° 


1.4 


"maximum" SO component 


MB3 


5 


10« 


1.4 


nt = - \ non-homogeneous SO comp. 


MB4 


5 


10' 


1.4 


nt = -2, "isothermal" SO comp. 


Tl 


5 


10" 


1.4 


fadi = 100, long switch-on 


T2 


5 


10' 


1.4 


temporary SO comp. (off at f = 100) 


Ql 


5 


10" 


1.3 


varied Toomre parameter 


Q2 


5 


10' 


1.5 




Q3 


5 


10" 


1.6 




Q4 


5 


10" 


1.7 




Q5 


5 


10" 


1.8 




Q6 


5 


10" 


1.9 












varied mass of the polar disk 


MDl 


5 


10" 


1.4 


6.0 • 10" Mq 


MD2 


5 


10" 


1.4 


6.0- 10" Mo,buta„,„ = 1.2 


MD3 


5 


10" 


1.4 


2.4 • 10'° Mq 


MD4 


5 


10" 


1.4 


3.7 • 10'° Mq, maximum disk 


Rl 


5 


10" 


1.4 


varied rotation profile, n=6 


R2 


5 


10" 


1.4 


n=3 


SI 


5 


10" 


1.4 


varied Eq. of state, y = 2.0 


S2 


5 


10" 


1.4 


7 = 1.0 (isothermal) 


S3 


5 


10" 


2.1 


7 = 1.0 (isothermal), Q^in = 2.1 


Gl 


5 


10" 


1.4 


larger grid (686x686) 


G2 


5 


10" 


1.4 


inner edge at R,r,in = 1 kpc 


G3 


5 


10" 


1.4 


as G2, but smooth transition of E 


G4 


5 


10" 


1.2 


as G3, but reduced Q„^\„ 



4. Models without an SO-like central component 

The aim of this paper is to study the stability of polar disks per- 
turbed by a central SO-like component. As a test of the code, 
but also for comparison with the models of the following sec- 
tions, we first describe in Sect. |4] the evolution of an unper- 
turbed gaseous disk (model A). In Sect.|5lwe perturb the polar 
disk by a central SO-like component. A list of the models de- 
scribed in this paper is given in Tabled 

In the first set of simulations, we investigate the evolution 
of an isolated galactic disk which has a central hole as de- 
scribed in the previous section. Such a disk is characterized 
by an exponential mass profile, a basically flat rotation curve 
and a Toomre parameter of g = 1 .4, a value slightly above the 
one required for stability against axisymmetric perturbations. 

Inspecting the Lagrange radii, i.e. radii containing a con- 
stant mass fraction, we see that the system is radially stable: the 
relative changes of the Lagrange radii are less than 0.01% over 
the whole integration time of f = 300 = 4.47 Gyr. On the other 
hand, the surface density profile shows the development of a 
pronounced two-armed trailing spiral structure until the end of 
the simulation (Fig.|3. The spirals are tightly wound in agree- 
ment with a short critical wavelength A^n of about 3 to 4 kpc. 
The most unstable wavelength is then about A„ ~ 1.5 -2 kpc 
for a gaseous disk (BT87) which reproduces roughly the wave- 
lengths visible in Fig.|3l The maximum amplitudes of the per- 
turbed density AI,(R, <p, t) = I.(R, ip, t) - I,(R, <p,t -0) are be- 




Fig. 3. Perturbation of the surface density profile at the end of 
the simulation of model A (f = 300). Contour lines are given 
in M0 pc~^ starting from 0. 1 to 0.9 M0 pc"^ with a diff'erence 
of 0.2 M0pc"^ between the contour lines. For comparison, 
the unperturbed surface density at a distance of 5 kpc is 37.2 
Mopc-2. 



low 1 Mqpc~ . Normalized to the initial densities, A2/S is 
much less than 1 %, which would not be detectable in observa- 
tions of real galaxies. 

Another way to visualize the nature and growth of insta- 
bilities in disks is to consider the azimuthal Fourier modes a,„ 
at different radii. They are calculated for m by (cf. also 
Laughlin et al. lT^ 



a„(R, t) 



27r Jo 



1(7?, (fi, t)e'""fdip 



Their phase angle is then given by 
5[a,„iR,t)] 



n(R, t) = tan 



5l[flm(^,f)] 



(14) 



(15) 



This phase angle is a measure for the spatial orientation of the 
related perturbation. E.g. a bar described by E(/?, ip) - I,oiR) ■ 
{1 +A cos[2{(p - a)]} has its major axis along (p - a 01 if = a+n, 
respectively. The corresponding phase angle for the m-mode is 
<pm(R) = I'na. Thus, the pattern speed for a mode m can be 
calculated by 



Q.,(R,t) = 



1 



„(R,t) 



dt 



(16) 



A measure for the global structure of a disk are the (radially 
integrated) global Fourier modes defined by 



CM = 



In 



a,n{R, t)RdR . 



(17) 



Mdisk is the mass of the disk within the radial range {Rio-x, ^hig]- 
The global Fourier amplitude is the modulus of C,„, whereas 
the global growth rate is given by the time derivative of the 
logarithmic Fourier amplitude, i.e. 



Jm = 



d{\n\Cjm 
dt 



(18) 
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Co 








X X X X 


^.1,1,1,1,1. 


A A a A ^ = 3 

□ □ □ □ 





50 100 150 200 250 300 

time 

Fig. 4. Temporal evolution of the global Fourier modes for 
m - 1 , . . . , 4 of model A. The mode m = describes the ra- 
dial matter redistribution. For details see text. 

It is also common to define a m = 0-mode which denotes the 
temporal evolution of the radial mass redistribution: 

Co = -^ f '"\m,t)-'^o(R)\RdR. (19) 

We adopted /?iow - 2 kpc and /?hig - 50 kpc, integrating over 
the entire disk. £(/?, f) and £()(/?) are the azimuthally averaged 
surface densities at time f, or at the start of the simulations, 
respectively. As long as the evolution is well described by a 
linear approach, no radial mass transport should occur. Thus, 
the increase of Co is a good indicator for the onset of non-linear 
effects. 

Fig- S shows the temporal evolution of the global Fourier 
modes for the "punched" disk without a central component. 
Starting from a low noise level, the even modes m - 2,4 
are the dominant ones describing a basically two-armed struc- 
ture. In agreement with predictions of linear stability analy- 
sis they grow exponentially on a similar, but long timescale of 

~ 20 ~ 300 Myr. In a late stage the m - 3-mode grows, 
too, whereas the m = 1-mode remains stable throughout the 
whole simulation. Until the end of the simulation the system 
never reaches the regime of non-linear growth or saturation. 
The amplitudes reached after 4.4 Gyr are so low that in real 
galaxies they would be practically invisible. 

The Cfl-amplitude shows a small initial radial mass redis- 
tribution which is caused by the small deviations from equi- 
librium due to the imprinted initial density perturbation. After 
quickly establishing a new equilibrium, the system remains ra- 
dially stable. 

The evolution of the phase angles gives more detailed in- 
formation on the structure formation. During the initial stage 
of the simulation no large scale pattern is discernible. However, 
after time f ~ 10 coherent features are detectable for the phase 
angle of the m = 2 mode which become more and more promi- 
nent indicating the weak two-armed structure of the system 
(Fig- El upper diagram). The structure grows from inside out, 
starting at a distance of about 4 kpc, well outside the inner edge 
of the computational grid. Therefore, it is very unlikely that the 
growing instability is caused by the sharp cut-off due to the 
central hole of the disk or by boundary effects. 
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Fig. 5. Evolution of the phase angle <p„ of model A for the m — 
2 mode during the initial stage (upper diagram) and the growing 
stage (lower diagram) (0°: white; 360°: black). 
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Fig. 6. Evolution of the radially averaged (R < lOkpc) pat- 
tern speed of the m - 2 mode (solid) and its mean deviation 
(dashed, lower curve) for model A. 

In a later stage when the linear growth is well established 
the tightly wound spiral structure is clearly visible in the phase 
angles (Fig.|5j lower diagram). These allow an accurate deter- 
mination of the intrinsic pattern speed of the two-armed spiral 
in the late stage as Qp,,„=2 ~ 1.07 a; 70kms"'kpc"' (Fig. |6}. 
^p.m=2 is the radial average of the local pattern speed Q.p{R, t) 
within the innermost 10 kpc. The small mean deviations of 
0.p{R, t) within this region indicate also the formation of a 
large-scale m - 2 pattern after t - 50. For this pattern speed. 
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Fig. 7. Schematic view of the geometry of the simulated polar 
disk. 



there are no inner Lindblad resonances and corotation is close 
to but inside the inner boundary of the grid. Similar to the m=2 
mode, regular patterns are also found for the m=3- and m-4- 
modes. Their pattern speeds are slightly smaller than that of 
the m=2 mode (by 4% and 7%, respectively). However, they 
are still too large to yield corresponding ILRs. 

5. Models with a central perturbation 

In this section we describe the evolution of "punched" disks 
perturbed by a central SO- or elliptical-like component as a 
model for polar disks. 

The perturbation due to the SO galaxy is given by the sta- 
tionary potential of a Ferrers' bar which has a density distribu- 
tion (BT87, Sect. 2.3) 



(20) 
(y/bf + 



9 / o\^h 

Ph(m ) = po ■ (1 - j ■ 0(1 - m) . 

is here the Heaviside function and = (x/a)^ 
(z/c)^ a modified radial coordinate. For most of our simulations 
we choose «/, - assuming a homogeneous mass distribution. 
In our reference model, the total mass of the central component 
is set to Mbar = 5-10' M0. 

The SO galaxy's major axis has a radial extent of 1 kpc, so 
the SO disk lies completely inside the central hole of the polar 
disk. The axis ratios are set to b/a - 1 and cja - 0.4 corre- 
sponding to an oblate spheroid. The minor axis points to the 
z-axis of the coordinate system. Thus, the "disk" of the central 
component is located in the jc-y-plane and the polar disk is in 
the y-z-plane. Hence, the disk plane of the central component 
corresponds to the abscissa of the polar disk snapshots, while 
the axis of the ordinate points to the z-axis, perpendicular to the 
disk plane of the central SO component (as shown in Fig.0. 

In order to avoid switch-on effects, the bar-like potential of 
the central SO component is turned on slowly on a timescale 
of fadi -2Q K 300 Myr. During that time the mass of the SO 
component is increased, while the same amount of mass is sub- 
tracted from the bulge/halo contribution. 

5.1. The reference model 

In order to start with a model resembling closely model A with- 
out a central component, we kept the rotation curve and the 
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Fig. 8. Rotation curve vv of a polar disk including a central 
component (model B). Shown are the contributions of the dif- 
ferent components (vhaio: halo, Vdisk: disk, vp: pressure and vso: 
SO component). 



-4 



-10 



1 1 1 1 1 1 1 1 1 








^Jj,xxxxxxxxx^x~^«f8^x X X □ 
X x^n A^^« 

X xDRD 


• • • • m=1 

X X X X 

" " " 'm = 3 
□ □ □ □ 










A A ^ ^ 




1 . 1 . 1 . 1 . 1 





50 100 150 200 250 300 
time 



Fig. 9. Temporal evolution of the global modes of model B for 
m - 1 , . . . , 4. The mode m - describes the radial matter 
redistribution. 



mass distribution in the polar disk constant, but added the cen- 
tral perturbation with Mbar = 5 10^ (model B). Due to the 
additional gravitational force generated by the central compo- 
nent, less mass is now in the halo. However, the velocity curve 
is still halo-dominated (Fig.|S}. 

The temporal evolution of the global modes (Fig. |5) dif- 
fers significantly from the simulation without a central pertur- 
bation. The perturbed disk develops quickly linearily growing 
even modes {m - 2 and m - 4), whereas the odd modes remain 
on their initial low level. After a time t ~ 120 a; 1.8 Gyr these 
modes become non-linear, which is reflected in an enhanced 
growth of the m = 0-mode. This radial mass redistribution is 
also shown by the Lagrange radii (Fig.llO>. The mass inside the 
10%-radius flows inwards, whereas the matter up to the half- 
mass radius expands. At the time t ~ 150 ~ 2.2 Gyr the m - 2- 
mode reaches its (first) saturation level of about 10%, while the 
m = 4-mode reaches a lower amplitude of ~ 10"^ - 10"^. At 
about the same time the odd modes begin to grow, but they are 
still linear at the end of the simulation at f = 300. Therefore, 
throughout the whole evolution the pattern looks mainly two- 
armed with a slight modulation introduced by the m - 4-mode 
(Figs.[TT]and[l2li. 
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Fig. 10. Temporal evolution of the Lagrange-radii containing 
0.1%, 1%, 5%, 10% (thick line), 20%, 50% (dashed), 75%, 
90% (thick) and 99% of the disk mass for model B. 
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Fig. 11. Perturbation of the surface density profile of model B 
during the linear regime (f = 60 = 890 Myr). Contour lines 
are given in M0pc"^ starting from 0.1 to 0.9 M0pc"^ with 
a difference of 0.2 M0pc"^ between the contour lines. The 
surface density at 10 kpc is about 10 M0 pc"^. Note the leading 
arms of the formed pattern. 



Fig. 13. Left, B-band image of NGC 4650A with spiral struc- 
ture indicated (Figure 10 of Arnaboldi et al. I1997> : right, the 
surface density distribution of Fig.^jseen in projection, as if 
the disk is tilted about the y-axis. 




6 
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Fig. 14. Evolution of the phase angles for the m -2 mode (0° 
white; 360°: black). 
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Fig. 12. Surface density profile in the saturation stage at the end 
of the simulation of model B (f = 300 = 4.47 Gyr). Contour 
lines are given in M0pc"^ starting from 10 to 100 M0pc"^ 
with a difference of 10 M0pc"^ between the contour lines. 
The maximum relative density perturbation AE/Zq is about 5.5 
in the knot at (0,-6)and 1 1.6 near the center 



It is remarkable, however, that the initial spiral pattern is 
leading (as Fig. shows) whereas the final spiral pattern is 
trailing like in the unperturbed model A. Moreover, the oc- 
curence of an initial leading pattern is found in all our models 
with a central SO component. It is especially independent from 
boundary conditions or other "technical" parameters. 

The temporal evolution of the surface density is depicted 
in Fig. [21 during the linear stage (until t - 120) a regu- 
lar tight leading spiral pattern emerges which becomes more 
and more pronounced. Later, large density variations show 
up which make the appearance temporarily ring-like. E.g. at 
f = 180 there is an almost closed ring at a radius of about 4 
kpc. Towards the end of the simulation a new quasi-equilibrium 
is established. This new mass distribution leads to more open 
trailing spiral arms in the saturation phase (cf. also Fig. I12> . 
When projected as if seen close to edge-on, the spiral qualita- 
tively resembles that seen in NGC 4650A (Fig. [Oj- We will 
show below that this open spiral is essentially independent of 
the central SO galaxy, which acts only as a 'seed' for its forma- 
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Fig. 15. Evolution of the logarithm of the surface density (in M0 pc"^) at different times (f = (upper left), 60 (upper right), 
120, 180, 240, 300 (lower right)) for model B. For a better visualization of the "punched" central area a cylinder with a constant 
level of 3.0 is drawn. 



tion. Obviously, the two-armed perturbation created by the cen- 
tral bar-like SO component induces a mainly two-armed struc- 
ture in the disk. This becomes clear in the diagram of the phase 
angles (Fig. ll4t . Immediately after the start of the simulation a 
regular, but stationary pattern is formed from inward out. This 
stationarity can be understood as the response of the flow to 



the stationary central perturbation. At a time f ~ 90 - 100 this 
large-scale tight leading pattern vanishes and a rapidly rotat- 
ing open trailing spiral is formed. The pattern speed (derived 
from the temporal change of the phase angles) at that time is 
about 0.9 ~ 60 kms"' kpc"' which is close to, but lower than 
the intrinsic pattern speed of the unperturbed disk. The tran- 
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Fig. 16. Temporal evolution of the amplitudes of the m - 0, 2- 
modes for the reference model B and a small central component 
(10'' M0, model MBl). 



Fig. 18. Temporal evolution of the amplitudes of the m - 0, 2- 
modes for the reference model B and a central component 
growing slower by a factor of 5 (model Tl). 




Fig. 17. Temporal evolution of the amplitudes of the m - 0, 2- 
modes for the reference model B and a maximum central com- 
ponent (10"' Mq, model MB2). 

sition from stationary pattern to rotating pattern occurs when 
the azimuthal force of the growing spiral structure begins to 
exceed the azimuthal force of the central component. They be- 
come comparable for the first time at about f ~ 85. Afterwards 
the azimuthal force of the spirals quickly outnumbers the con- 
tribution of the bar-like central component, making the central 
perturbation unimportant for the further growth of the spiral 
structure. 

5.2. Variation of the central component 

In order to test the influence of the central component we varied 
the mass of the central component (models MBl and MB2), its 
switch-on characteristics (model Tl) as well as the duration of 
the central perturbation (model T2). 

Mass of central component. Simulations with a bar-Uke 
SO component that is five times less massive (1 ■ 10^ M0: 
model MBl) and with a "maximum SO component" of 
10'° M0 (model MB2) have been performed. The smaller 
mass leads to a delayed growth of the spiral modes of the disk 
by about f ~ 25 * 370 Myr (Fig.[^. In case of the more mas- 



Fig. 19. Temporal evolution of the amplitudes of the m - 0, 2- 
modes for the reference model B and model T2 in which the 
central component is switched off between f = 80 and f = 100 
(1.2-1.5 Gyr). 

sive central pertubation there is again a small temporal shift in 
the behaviour of the amplitudes during the linear stage (Fig. 
I17> . The amplitudes run now ahead by about 100-150 Myr In 
both cases the growth rates and the final saturation levels are 
unaffected by the varied mass of the central component. 

We also varied the density profile of the central component 
by choosing = (MB3) and nt - -2 (MB4). However, 
even the steeper profile did not change the numerical results 
significantly. 

Timescale of perturbation due to central component. 

Another test of the influence of the central component is pro- 
vided by a different temporal behaviour of the perturbing bar- 
like potential. In model Tl the switch-on time fadi was increased 
by a factor of 5 to fadi - 100. Similar to the model MBl with a 
less massive central perturbation the growth of the main ampli- 
tudes is now delayed by about 400 Myr ('Fig.ll8>. The growth 
rates and the saturation levels, however, remain unaffected. 

Just as in the previous models, even a complete removal of 
the SO component does not alter the evolution strongly. Figll9l 
shows that the disk quickly evolves on its standard evolution- 
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Fig. 20. Temporal evolution of the amplitudes of the m - 2- 
modes for the reference model B (1.2- 10'° M0, plus) and mod- 
els with different disk masses: MD2 (6.0- 10^ M0, triangle) and 
MD3 (2.4 ■ 10'° Mq, filled circle). 




50 100 150 200 250 300 

time 

Fig. 21. Temporal evolution of the amplitudes of the m - 
2-modes for the reference model B and models Ql 
. . . Q6 with different initial minimum Toomre parameters 
(Q=l. 3,1.6,..., 1.9). 



ary path ending up in the same saturation stage, even if the SO 
component is switched off before the density perturbations in 
the disk reach the non-linear regime. The perturbation induced 
by the SO component does not disappear when the aspheric- 
ity of the central perturbation is switched off, and the modes 
continue to grow. When the gaseous component is dynamically 
substantially hotter, we find that the excited modes are erased 
after the perturbation is switched off. 

From these calculations we conclude that the modal growth 
in the reference model B decouples from the initial central per- 
turbation after time t ~ 100. The central component acts as 
a seed for the growing modes, but once these modes are ex- 
cited, they decouple from the seed growing on a timescale of 
1-2 Gyrs. Also other perturbations, e.g. due to gas accretion, 
might act as such a seed. 

5.3. Mass of the polar disk 

The influence of the mass of the polar disk is studied by a small 
series of simulations starting from half of the mass of the refer- 
ence model (MDl, MD2) up to the maximum polar disk mass 
of 3.7 ■ 10'" M0 (MD4). A disk with half the mass of the ref- 
erence polar disk (or about 1/6 of the maximum polar disk) is 
very stable (MDl): after an initial adjustment the amplitudes 
of all modes remain almost constant on a very low level. The 
dominant m - 2 and m = 4 modes are still deep in the linear 
regime. Only a weak and very tightly wound leading spiral is 
formed similar to that found in the early stage of the reference 
model. Model MD2 has the same lowered disk mass as model 
MDl, but a reduced minimum Toomre parameter Qmin - 1-2. 
Still the disk is rather stable (Fig. I20> and the dominant even 
modes have the same amplitude as in MDl. The main differ- 
ence to MDl is that the odd modes rise very slowly, but after 
4.5 Gyr their amplitudes are still well below those of the even 
modes. 

Increasing the disk mass to twice the mass of the reference 
model, destabilizes the polar disk (MD3, Fig.l20>. The growth 
rates of them - 2 and m - 4 modes are doubled compared to 



the reference model. The growth rates of the m - 1 and m - 3 
modes are very similar to those of the even modes. The overall 
appearance within the first 2 Gyr, however, is dominated by 
the even modes because their early amplitude is raised by the 
perturbation due to the central SO component. Increasing the 
mass up to the maximum polar disk (MD4) gives a result very 
similar to the model with doubled mass MD3. 

5.4. The Toomre parameter and the rotation curve 

An important parameter for the stability of disks is the Toomre 
parameter. In a series of models we varied the initial minimum 
Q-value between 1.3 and 1.9. Our reference model B was char- 
acterized by 2min = 1.4. Fig. 12 II demonstrates the strong de- 
pendence of the growth rate on Q during the linear stage: hot- 
ter disks are more stable. Close to the reference model (e.g. 
comparing models with g = 1.3 and Q = 1.5) the growth rate 
varies less strongly than in models with larger Q-values (e.g. 
the range Q = 1.5 to Q = 1.7). Moreover, the disks show a 
transition to complete stability when the Toomre parameter ex- 
ceeds a value between 1.7 and 1.8. This is in good agreement 
with analytical predictions for the transition to stability to be at 
Qc ~ 1.73 (in case of flat rotation curves, see e.g. GBOO). The 
saturation levels of the dominant m - 2-modes are - in case 
of unstable configurations - rather similar, i.e. of the order of 
about 10%. 

In the models Rl and R2 we tested the influence of a 
broader transition region between the rigid rotation domain and 
the flat part of the rotation curve by setting n - 6 (Rl) and 
n = 3 (R2). Model Rl develops qualitatively similar to the 
reference model B. However, the growth rates of the modes 
are 30% lower and, therefore, the non-linear regime is first 
reached after about 3.5 Gyr. In model R2 the transition region 
is even broader. Still a growing instability is found, but it grows 
on a longer timescale. The structures dominated by the lower 
even modes become non-linear at about t x 300 ~ 4.5 Gyr. 
Following the linear stability analysis of Polyachenko et al. 
(i 1997 J the critical Toomre parameter Qc (for which the system 
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Fig. 22. Temporal evolution of the m=2-amplitude for different 
equations of state: y=5/3 (reference model B, plus sign), 7=2 
(model SI, open box) and 7=1 with 2min=l-4 (model S2, open 
triangle) and 2min=2.1 (model S3, filled circle). 
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Fig. 23. Temporal evolution of the amplitudes of the m - 0, 2- 
modes for the reference model B (270x270) and a model Gl 
with increased grid size (686x686). 



is stable against all perturbations; see their Eq. (36)) drops from 
Q^= V3 ~ 1.73 for flat rotation curves to = V372 * 1.22 
for the rigid rotation region. Thus, a broader transition region 
results in a reduced beyond the transition radius /?flat com- 
pared to the reference model, and the system is expected to be 
more stable, in agreement with our simulations. 

5.5. Equation of state 

In the models SI to S3 we changed the equation of state by 
varying the adiabatic exponent y in the polytropic equation of 
state. In model Sly was set to 2 which is often adopted to 
mimic stellar disks (e.g. Orlova et al. '2002V Still such a disk 
is not stable, however the growth of the Fourier amplitudes is 
delayed and smaller ("Fig. l22t . Such a stabilizing behaviour can 
be understood from the stiffer equation of state for y = 2 which 
hinders perturbations from growing. 

Correspondingly, the disk is expected to be more unstable 
when the equation of state is less stiff as e.g. in the isothermal 
case y = 1 . Here cooling is taken to be so efficient that all the 
heat of compression is dissipated immediately, so that the disk 
is maintained in an isothermal state, as suggested by Lodato 
& Rice (2005). Accordingly, we see that model S2 reaches the 
non-linear regime in about half the time of the reference model. 
However, the destabilizing effect of the reduced stiffness can 
be partly compensated for by an initially hotter disk, as model 
S3 with a minimum Toomre parameter of Qmi,, =2.1 demon- 
strates. The cases y = 1,5/3 and 2 are expected to bracket the 
expected cooling behaviour in real disks. 

5.6. Miscellaneous 

We tested different variations of "technical" parameters which 
did not result in substantial quantitative deviations from the 
reference models. E.g. we varied the numerical timestep cri- 
terion, we compared with models using different amounts of 
artificial viscosity or we changed the grid size. As an example 
Fig. l23l demonstrates the insensitivity of the numerical results 



on the grid size. Though the resolution was more than doubled, 
the amplitudes characterizing the dominant mode or the radial 
mass redistribution evolve qualitatively in the same way and 
quantitatively almost identical. 

In order to test the influence of the inner boundary we rerun 
the unperturbed model A, but with an inner edge at R^in - Ikpc 
(models G2 and G3). In both models we kept the total disk mass 
constant. In model G2 the exponential disk was extended down 
to a radius of 1 kpc, while in model G3 the surface density was 
exactly exponential outside a radius of 2 kpc (like in model A), 
but had a smooth transition to zero from a radius of 2 kpc down 
to the inner edge. Both models were more stable than model A. 
In model G2 the growth rate of the dominant m = 1-mode is a 
factor two smaller than the growth rate of the dominant mode 
OT = 2 in model A. Model G3 shows no significant growth at 
all. In model G4 the same smooth surface density distribution 
as in model G3 was used, but the minimum Toomre parameter 
was reduced to 1.2. In that case, the model becomes violently 
unstable after 4 Gyr. This means that the stabilizing effect of 
a smoother mass distribution can be compensated by a slightly 
reduced minimum Toomre parameter. Taking into account that 
the velocity dispersions within the polar disk are poorly known 
(if at all), the uncertainties related to the mass distribution at 
the inner edge of the disk are less critical, unless the profile of 
the Toomre parameter can be constrained in more detail. 

6. Discussion and summary 

We investigated the stability properties of polar disks by 2- 
dimensional hydrodynamical simulations of flat disks per- 
turbed by a non-rotating bar-like potential mimicking a central 
SO-like component. The properties of the polar disk and the 
central perturbation were derived from observations of the pro- 
totypical polar ring galaxy NGC4650A. We modelled the polar 
ring by a disk with a "punched" central hole and an exponential 
mass distribution embedded in a halo with a mainly flat rota- 
tion curve. We modelled the polar gas using a polytropic equa- 
tion of state (EOS), and did not consider a multi-phase ISM or 
energy feedback from star formation in our calculations. Both 
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will probably influence especially the late non-linear stages of 
the evolution, when cool clumps can form and star formation 
sets in. E.g. star formation is expected to stabilize the disk due 
to its energy feedback to the ISM, but a population of cool 
dense clouds is destabilizing (e.g. Jog & Solomonil984 Shen 
& Lou2Q04). However, the timescale for the disk to become 
unstable should not be affected strongly. In order to study dif- 
ferent kind of gas behaviour we varied the adiabatic exponent 
covering the isothermal case (y = 1) up to the very stiff "stel- 
lar" case (y - 2). Though the numerical results depend on the 
actual EOS, its influence also strongly depends on the Toomre 
parameter of the disk, in the sense that a hotter disk can com- 
pensate for a less stiff EOS. 

Without the central perturbation, the disk is fairly stable 
over the whole simulation time of 4.5 Gyr. Two-armed and 
four-armed modes grow, but their growth rates are rather small. 
The corresponding Fourier amplitudes of 10 are even after 
4.5 Gyr deeply in the linear regime and, therefore, basically 
undetectable to observation. 

Exciting the perturbation by a central bar-like component 
representing the SO galaxy changes the evolution of the polar 
disk drastically. Though the perturbation was switched on adi- 
abatically, the Fourier amplitudes reach a level of 10~^ more or 
less as soon as the SO potential is switched on. Tightly wound 
leading spirals are formed; though their flow pattern remains 
stationary, their amplitude is growing. In our simulation the 
azimuthal component of the gravitational force exerted by the 
spiral perturbations becomes comparable to that of the central 
SO at about 1 Gyr. The evolution of the polar disk then begins 
to decouple from the central component. At about 1.8 Gyr a 
substantial radial mass redistribution sets in, indicating that the 
evolution becomes non-Unear. The tightly wound spiral disap- 
pears and, finally, a much more open, but rotating, two-armed 
trailing spiral develops. 

The general behaviour turned out to be rather robust against 
variations of the central component (e.g. its mass or the tem- 
poral behaviour of the perturbation). Switching off the non- 
axisymmetric component of the imposed external potential, 
even if it is done well before the nonlinearity has developed, 
has no apparent effect on the state at the end of the simulation. 

The most important parameter for the evolution of the po- 
lar disk is the Toomre parameter In case of a high Q exceeding 
1.7, only tightly wound spirals form, driven by the barlike per- 
turbation of the SO disk. Such disks are too hot to be unstable 
against the induced spiral perturbations, and the tight spirals 
disappear when the central component is switched off. This is 
exactly the reason usually given for the absence of spiral struc- 
ture in the stellar disks of SO galaxies: the disk is too hot to re- 
spond coherently and support a self-excited global spiral mode 
(e.g. Chapter 18 of GBOO). In case of lower Q (1.6 or less) the 
central component just acts as a seed of instability. In the early 
stage a tightly wound leading spiral is formed like in the refer- 
ence model. This has a strong resemblance to the spiral in the 
inner polar disk of NGC 2787 (see e.g. Erwin & Sparke 2D03 1. 

Later on, the spiral structure decouples and becomes an 
open trailing spiral during the saturation stage. The timescale 
for reaching the saturation stage slightly depends on Q ranging 
from 2.2 Gyr (Q = 1.3) to 3 Gyr (Q = 1.6). Taking a Fourier 



amplitude for the m - 2-mode of 1 % as the margin for a de- 
tectable spiral structure, the timescale for our reference model 
of a cold polar disk to become unstable is about 1-2 Gyr The 
transition to the saturation stage takes then another Gyr These 
timescales are remarkably independent of the mass of the cen- 
tral component, and depends mainly on the mass in the polar 
disk and the Toomre parameter If the mass is halved, the insta- 
bility disappears unless the Toomre parameter is lowered, and 
the trailing spiral saturates only after 4.5 Gyr. If the polar-disk 
mass is doubled, growth rates are also approximately doubled, 
but the saturated spiral looks similar Again, studies of stellar 
systems have long noted that relatively cool disks containing 
a large fraction of the system's mass are likely to be unstable 
to forming a spiral pattern: see e.g Binney & Tremaine J1987> 
for a summary. The spirals formed in that stage are all trailing. 
The open spiral in NGC 4650A described by Arnaboldi et al. 
( 1997) and Gallagher et al. (2002i might be an example of a 
late stage of such an evolution. 

On the basis of these computations, we expect an initially 
axisymmetric polar disk resembling that in NGC4650A to be 
stable at least over 1-2 Gyr, or longer if it is hot. When they 
become unstable, they first form tightly wound leading spirals, 
before they rearrange their mass building up an open trailing 
spiral. Hence, the variety of observed structures in polar disks 
might imply a large range of ages. 

In order to predict the long-term dynamical evolution of a 
polar disk, it is essential to determine its Toomre parameter. 
Hence, in addition to the rotation curve and the mass distribu- 
tion, the velocity dispersions must be known. 

Additionally, it would be very interesting to test whether 
observed polar disks generally have some substructure and if 
so, how many of them show tightly wound or open spirals. 
Moreover, it would be very important to check whether there 
are tightly wound spirals which are indeed leading. Such a mea- 
surement, however, is very difficult due to the small Fourier 
amplitudes in the stage when leading arms are expected. In the 
transition stage when the leading spirals become trailing, the 
arms are still tightly wound. Hence tightly bound arms can be 
trailing or leading depending on their evolutionary stage. 

Another major problem is that polar disks are most easily 
visible when they are almost edge-on. However, then it is dif- 
ficult to analyse the substructure in the disk (see e.g. Figs. 3a- 
3d in Whitmore et al.'s ( 199()_) catalog of polar-ring galaxies). 
However, deep images of a selected sample of inclined polar 
disks might allow for such a measurement. 

Bekki & Freeman (2002.) have given an alternative expla- 
nation for the spiral in the polar disk of NGC 4650A. The per- 
turbation by a triaxial dark matter halo tumbling slowly about 
an axis perpendicular to the polar ring can cause an open spiral 
structure well outside the region where we expect the polar disk 
to be self-gravitating. They emphasize that self-gravity is less 
important for their structure formation mechanism, because 
they produce mainly two-armed kinematic spirals. However, 
the rotation of the dark halo must be finely tuned to produce 
the observed open spiral. 

We have demonstrated that the (unavoidable) perturbation 
given by the central (SO) component of a polar ring galaxy 
embedded in a spherical dark matter halo component can be 
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sufficient to create open trailing spiral arms well outside the 
central SO galaxy. The self-gravity of our polar disk is essential 
for creating pronounced spiral structures with amplitudes in the 
non-Unear regime, which need not necessarily to be dominated 
by even modes. On the other hand, in low mass polar disks it 
is very difficult to form outer spirals by a central perturbation. 
We may see a tight trailing spiral forced by the perturbation of 
the central oblate galaxy, or no spiral pattern at all. 
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